Precursors of order in aggregates of patchy particles 
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Abstract 

We study computationally the local structure of aggregated systems of patchy particles. By 
calculating the probability distribution functions of various rotational invariants we can identify 
the precursors of orientation order in amorphous phase. Surprisingly the strongest signature of 
local order is observed for 4-patch particles with tetrahedral symmetry not for 6-patch particles 
with the cubic one. This trend is exactly opposite to their known ability to crystallize. We relate 
this anomaly to the observation that a generic aggregate of patchy systems has coordination number 
close to 4. Our results also suggest a significant correlation between rotational order in the studied 
liquids with the corresponding crystalline phases, making this approach potentially useful for a 
broader range of patchy systems. 
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The field of colloidal and nanoparticle self-assembly has dramatically changed over the 
past decade due to introduction of novel classes of particles and interactions between them. 
This includes highly selective DNA-mediated interactions PQ-[3], use of wide variety of parti- 
cle shapes and so-called "patchy" colloids [I]- [6], as well as combinations of these approaches. 
The patchy particles have patterns of chemically distinct regions on their surfaces that results 
in directional (covalent-like) interactions. This opens an appealing prospect of "program- 
ming" the symmetry of a desired structure with the symmetry of the particle. For instance, 
colloids with tetrahedral arrangement of patches have been widely studied theoretically [H]- 
[TT] . as a potential platform for self-assembly of diamond lattice, one of the best candidates 
for photonic band gap materials [7j. These studies were in part motivated by experimental 
demonstration of patchy colloids with tetrahedral and other symmetries [5j. Most recently, 
these experimental techniques evolved even further due to functionalization of patches with 
DNA and resulting selectivity of interactions [6]. 

While the equilibrium phase diagrams of such patchy colloids have been extensively 
studied computationally, the experimental self-assembly of crystalline morphologies will 
unavoidably be limited by slow kinetics. By analogy with self-assembly of colloids and 
nanoparticles isotropically functionalized with DNA, one might expect the system to form 
a random aggregate initially, and possibly be transformed to a crystal upon annealing. In 
this paper we focus on the relatively early stage of self assembly, and analyze the precursors 
of crystallinity in a (mostly) random aggregate of patchy particles. We do this by employ- 
ing the set of bond order parameters [15J which have been successfully applied to a great 
variety of physics problems. The parameters are commonly used to determine the degree of 
crystallinity of a system, as well as the crystal morphology. In our case, we apply it to find 
traces of ordering in liquid phase, which may result in crystal formation upon annealing. 
This kind of analysis can potentially determine the morphology which is preferred kinet- 
ically, rather than energetically. For instance, it is well known that the space of possible 
structures for patchy particles with tetrahedral symmetry is highly degenerate. At least two 
diamond morphologies, cubic, and hexagonal have nearly the same free energy. As a result, 
in majorities of studies (with a rare but noteworthy exception [11] ) neither of the crystals 
form spontaneously. In such a case, the structure in actual experiments may be selected 
kinetically. 

Several models have been used to simulate patchy colloids. Most common are anisotropic 
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Lennard- Jones potential [10], and Kern-Frenkel model [12] , [8] , [9] . In those models, the 
patch geometry can be tuned independently of the interaction range, and the degree of 
directionality is an important parameter of the system. Since the dependence of the system 
behavior on the patch size is not a focus of our study, here we use a simpler model with point- 
like patches interacting via Gaussian potential. The range of the potential automatically 
determines the degree of directionality of the bond. Physically, this is a reasonable model to 
describe nanoparticles with locally grafted DNA molecules, in which case both patch size 
and range of interactions is determined by DNA gyration radius. 

Numerical method Our numerical algorithm realize the description of a patchy particle as 
a rigid body. We represent the motion of a patchy particle of the radius R as a combination of 
the displacement of its center and rotation around some axis, passing through its center. We 
study two systems with 4 patches (4pch) and 6 patches (6pch). At the initial time moment 
orientation locations of four patches a^(0), k = 1, 2, 3, 4 for 4pch particles with tetrahedron 
symmetry with respect the center of the j-th particle are sl^' 2 \o) = f±<y/|ii!, 0, y^R), 

sl^' A \o) = ^0, ±y|i2, — yJ^Rj. For 6pch system with cubical symmetry at initial time 
moment patches are located at points a.j\o) = (±R, 0, 0), a^ 3,4 ^(0) = (0, ±R, 0), a.f' 6 \o) = 
(0,0,±R). 

The displacement of the center of the j-th particle at the time moment t is described 
by the vector Tj(t). The orientation of j-th particle at the time moment t is given by the 
unit quaternion Aj(t). That quaternion may be represented in the form (see, e.g., [13J) 
Aj(t) = [cos(0j(t)/2),sin(0j(t)/2)rij(t)], where the unit length vector |nj(t)| = 1 describes 
the axis, passing through the center of the particle and <frj is the angle of rotation around this 
axis. We also introduce conjugated quaternion Aj(t) = [co8(<f>j(t)/2), — sm(<f>j(t)/2)nj(t)]. 
Finally, the location of the k-th patch of the j-th particle at the time moment t is given 
by formula a.j(t) = Tj(t) + Aj(t) eg) a^(0) £g> Aj(t) where ® denotes the quaternion's prod- 
uct. In our model cores of patchy particles repel each other with standard Lennard- Jones 
potential smoothly truncated at the distance 2R [14J with the interaction distance a = 2R 
and interaction strength e = 1. Patches of different particles attract each other with the 



Gaussian potential Ua^aq!! 1 ') = U p exp 



-(af) 2 /2W 2 



where a^ = gq 1 ' — a-- is a vector 
connecting a patch I of particle j and a patch k of particle i, W = 0.2 is the half-width 
of the interaction and U p is the strength of the interaction. Knowing the set of all dis- 
placements and orientations of particles {rj, Aj} we can compute a set of total forces and 
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torques {Fj, Mj} acting on every particles. We solve the following set of kinematic equations 
numerically using the velocity Verlet algorithm 

irj(t) = Fjdrj, Aj})/m ( rj(t) = v,(t) 

= M,({r,, A,})// ' \ Aj(t) = i Wj (f) <g> A&) ' 

where Vj is the velocity of the j-th particle and Uj is its angular velocity, Vj and Uj are linear 
and angular accelerations, respectively. We use normalized units: the radius of particles R = 
1, the mass is m = 1, the moment of inertia of a solid sphere / = |mi? 2 = 0.4. To simulate 
the interaction with solvent we add Langevin terms — 7V J -(t) + £j(t) and —^R 2 ojj(t) + 
Cj(t) to forces and torques, respectively, where 7 = 6nr]R is the friction coefficient for 
solvent viscosity rj, £j{t) and £j(t) are thermal noise terms with delta-correlated components 
(ti(t)Zj(fj) = 2 1 k B T5 lJ 5 a ^5 t>tl , (Q'{t)Cf(f)) = liR 2 k B T5 l>j 5 a ^5 t , tl . In our simulations, 
ksT = 1 and 7 = 10, therefore for times t ^> I/7 = 0.1 the dynamics of a particle is 
Brownian. We simulate the system of iV = 10 3 spherical particles of radius R = 1 in 
a cubic cell of size L = 48 with periodic boundary conditions. The volume fraction is 
= 4/3nR 3 N/L 3 ~ 0.04. We can tune the phase state of this system by varying the 
strength U p of interaction potential between patches. 

Structural properties To define the local structural properties of the system we use the 
bond order parameter method [15J, which has been widely used in the context of condensed 
matter physics |15j, hard sphere systems [HUE]) complex plasmas [T8H20] . colloidal sus- 
pensions [21], granular media etc. Within this method the rotational invariants of rank I of 
both second qi(i) and third Wi(i) order are calculated for each particle i in the system from 
the vectors (bonds) connecting its center with the centers of its N nn (i) nearest neighboring 
particles 
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where qi m (i) = N nn (i) 1 Y^fZi^ ^/m( r ij), Yim are the spherical harmonics and = r.j — r^- 
are vectors connecting centers of particles i and j. We note, that the bond order parameters 
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wi oc qf, so, in general, these parameters are much more sensitive to the local orientational 
order in comparison with qi. Here, to define the structural properties of the patched particles, 
we calculate the rotational invariants qi, wi for each particle using the fixed number iV nn of 
the nearest neighbors (NN): iV nn = 4, 6 for the first shell of 4pch and 6pch systems; second 
shell of both systems has N nn = 12. The first shells of the ideal 4pch and 6pch systems 
have cubic diamond (CD) and simple cubic (SC) lattices, respectively. The second shell has 
face centered cubic (FCC) lattice for both types of ideal patchy systems. We note, that 
formation of hexagonal diamond (HD) with hexagonal close packing (HCP) second shell is 
also possible for the patched system, at least kinetically. 

The values of the different rotational invariants q\ and w\ for the perfect patchy crystals 
(for both first and second shells) are shown in Table 1. A particle whose coordinates in 
the 4-dimensional space (q^, q$, W4, Wq) are sufficiently close to those of the ideal lattice is 
counted as solid-like particle. By calculating the bond order parameters for the second shell 
it is easy to identify the disordered (liquid-like) phase as well. 



TABLE I: Rotational invariants for the perfect patchy crystals 



system structure 


94 96 


W4 


w 6 


4pch CD, HD (1st shell, 4NN) 


0.509 0.628 


-0.159 


-0.013 


6pch SC (1st shell, 6 NN) 


0.76 0.35 


0.159 


0.013 


4,6pch CD, SC (2nd shell, FCC) 


0.191 0.575 


-0.159 


-0.013 


4pch HD (2nd shell, HCP) 


0.097 0.485 


0.134 


-0.012 



Figure [T] and [2] show the probability distribution functions (PDFs) of bond order pa- 
rameters g 4 (a,b) and w 4 (c) at different strengths U p of the interaction potential for the 
first (a) and second (b,c) shells. Figure [I] and [2] correspond to the systems with 4 and 6 
patches, respectively. The MD simulations cover both uncorrelated (gas-like) and strongly 
coupled (liquid-like) phases. The PDFs show how the structural properties of the ensemble 
of patched particles vary with increase of interaction strength U p . At low U p (blue lines) 
both 4 and 6 patched systems are completely dispersed, and the PDF plots correspond to 
the isotropic distribution (of uncorrelated system). Increase of U p results in formation of 
aggregates of patchy particles; in the final state, nearly all particles consolidate into a few 
big clusters. 
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FIG. 1: (Color online) Patchy system 4pch. Probability distribution functions (PDFs) versus 
(a,b) and W4 (c) at different strengths U p of the interaction potential for the first (a) and second 
(b,c) shells. Cumulative distributions of the PDFs are also plotted to quantify the liquid-solid 
transition in the system. The curves are color-coded by U p value. Inset (a) shows first (green) and 
second (red color) shells of the perfect crystalline particle. Insets (b) show initial gaseous (b, left) 
and final gel-like (b, right) particle distribution over space. Particles are color-coded by q§ value. 

Review of the plots for the first shell reveal a remarkable and counterintuitive result. 
There is a strong evidence of local tetrahedral ordering for 4pch particles, represented by a 
spike at = 0.5, but very weak order for particles with cubic symmetry. This goes exactly 
contrary to the known crystallization properties of these systems: as we have discussed 
earlier, 4pch particles are notoriously hard to crystallize, as opposed to 6pch. The second 
shell PDF of shows a pronounced but wide peak centered around = 0.2, which is very 
close to the value expected for ideal FCC lattice. Note that the second shell of both SC 
and CD lattices has FCC symmetry. It is therefore possible that this is peak is a precursor 
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FIG. 2: (Color online). Patchy system 6pch. Probability distribution functions (PDFs) versus 
(a,b) and W4 (c) at different strengths U p of the interaction potential for the first (a) and second 
(b,c) shells. Cumulative distributions of the PDFs are also plotted to quantify the liquid-solid 
transition in the system. The curves are color-coded by U p value. Inset (a) shows first (green) and 
second (red color) shells of the perfect crystalline particle. Insets (b) show initial gaseous (b, left) 
and final gel-like (b, right) particle distribution over space. Particles are color-coded by q§ value. 

of the future crystalline order. In order to confirm this correlation, additional studies of a 
broader class of systems will be needed. If that is the case, it will imply that 4pch system 
prefers CD structure over HD, at least kinetically. Note that W4 does not show any clear 
signature of FCC order, but that may be due to the fact that it is a higher-order invariant 
than g 4 . In general both and W4 PDFs look remarkably similar for both systems, which 
may originate from the similarity between the second shell structure of SC and CD, but may 
also be a generic property of amorphous aggregates of patchy particles. It noteworthy that 
6pch exhibits very sharp but weak peak at W4 ~ 0.13. This value is close to the one for HCP 
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FIG. 3: (Color online) (a), (b) The order parameters, characterizing the crystallization of the 
system with 4pch (a,c) and 6pch patches (b,c) and associated with the cumulative PDFs versus q$ 
value are plotted as a function of the particle interaction potential strength U p for both first (FS, 
open triangles and squares) and second (SS, filled triangles and squares) shells. Panel (c) shows 
the normalized maximal cluster size S max /N (solid squares and triangles) and average coordination 
number Z (open squares and triangles) as functions of U p . Solid green lines correspond to analytical 
fitting discussed in text. 

order, but its actual origin is unclear at the moment. The HCP order is inconsistent with 
the cubic arrangement of the patches, and W4 value alone does not allow for an unambiguous 
interpretation. 

The cumulative distributions C q {q4) and C w {w±) are also presented in Figures 
One can use the half-height positions for these curves as order parameters that characterize 
the local orientational order. For instance, Q4 is denned based on cumulative distribution 
of C g (g 4 ) as J"^ 4 C q {q^)dq^ = 1/2. Its value for both first and second shells are plotted in 
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Figure [3] as a function of U p for both patchy systems. Note that the completely isotropic 
distribution corresponds to non-zero values of due to a finite number of particles in each 
shell. The deviation from that values characterizes the degree of orientational correlations 
in the system. We observe a clear crossover from isotropic gas phase to strongly correlated 
liquid at high values of U p . 

In order to understand the reason for the striking difference in the first-shell orientational 
order between the two systems, we will now consider their topological properties. Two 
patches are defined to be bound if they are closer than 2W from each other (where W = 0.2 
is the width parameter of the Gaussian potential). This allows us to construct clusters of 
connected particles and compute the fraction of particles that belong to the maximum clus- 
ter, S max /N. In addition, we determine the mean coordination number Z, i.e. the average 
number of neighbors to which a particle is connected. Both properties are presented in Fig- 
ure [3] together with cumulative bond orientation parameter Q% discussed above. Naturally, 
the signature of aggregation appears simultaneously on all the plots. Remarkably, the co- 
ordination number in both systems exhibits saturation at value Z « 4, despite the particle 
having dramatically different design and number of patches. This observation is the key to 
understanding the difference in first shell ordering. Indeed, this value of Z means that 4pch 
particles in random aggregate have almost the same number of bonds as in ideal diamond 
crystal, either cubic or hexagonal. Hence, the first shell exhibits strong tetrahedral ordering. 
On the other hand, at Z = 4 the 6pch system is far from connectivity of the corresponding 
simple cubic lattice, and the orientational order is far less pronounced. 

The dependence of the largest cluster size on U p can be also related to the value of 
Z. The concentration of free particles depends exponentially on the chemical potential of 
those belonging to the aggregate, fi = ZU p /2 + const. The size of the biggest cluster can 
be therefore estimated by subtracting the gas-like fraction of the system from the total: 
Smax/N — 1 — ex P ( — Z (Up — Uq) /2). This formula with Z=4 well describes the numerical 
results, as shown in Fig. 3(c). The relative shift of the two plots can be attributed to an 
additional entropy of the 6pch liquid. This entropy, AS ~ 4/cb per particle reflects a much 
larger number of ways in which 6-patch system can be arranged into into Z = 4 network. 

What determines the coordination number of the amorphous aggregate? In the limit of 
infinitely short range of interactions between patches, and hard-core interparticle repulsion, 
the bond orientation and relative positions of the two bound particles would be completely 
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restricted. The only remaining degree of freedom would be rotation about the direction of the 
bond. This means that each perfect bond freezes 5 translational and orientational degrees 
of freedom (or 2d — 1 in d dimensions). Since the total number of degrees of freedom for 
N particles is d(d + 1) N/2, the system becomes completely rigid and incapable of creating 
new bonds for coordination number Z rigi d = d (d + 1) / (2d — 1) = 2.4 for d = 3. This result 
represents conceptually important but somewhat unrealistic case from both experimental 
and computational points of view. Both interactions between patches and interparticle 
repulsions have finite range which means that bonds are not infinitely rigid. Interestingly, out 
of 5 degrees of freedom that a perfectly rigid bond would suppress, not all are equally affected 
in the case of finite interaction range. If Ai is the typical range of intra-patch attraction, 
and A2 is that of interparticle repulsion, a single bond results in an angular confinement of 
a particle within solid angle 5Q ~ 4ti (Ai + A2) / R. This means that each of the angular 
coordinates is constrained much weaker in relative terms, than the translational degrees of 
freedom of the bound patches, 59 ~ ^/ ( A i+ A z) ^> ^1 jf we now re p ea t the above counting 
argument by only assuming each bond to suppress d translational degrees of freedom, we 
obtain a new estimate of the coordination number of the amorphous aggregate Z* = d+1 = 4 
for d = 3. This value of Z is indeed consistent with our results, and also plays a prominent 
role in a other important problems. For instance, it corresponds to isostatic packing of 
objects with infinite friction coefficient. 

To summarize, the random aggregation naturally results in a liquid with coordination 
number close to Z* = 4. This is close to maximum connectivity of 4pch particles, but is 
substantially below that for 6pch system. As a result, the bond orientational order is much 
more pronounced in the first shell of particles with tetrahedral symmetry than for those with 
cubic ones. Paradoxically, this also means much lower driving force towards crystallization in 
the tetrahedral case. This complements the well known problem of degeneracy of the ground 
state of the system: not only is there a competition between energetically very similar cubic 
and hexagonal diamond, but both of them have very little advantage in connectivity over 
generic aggregate with Z rs Z* = 4. The 6pch system on the other hand has a clear energetic 
incentive to form a well coordinated SC crystal. Altogether this explains both why the 
diamond is so hard to self-assemble and also the anomalous orientational order that we report 
for 4pch liquid. In addition, our analysis of the second shell organization reveal signatures 
consistent with FCC ordering in both system, which is indeed expected in the corresponding 
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crystals, SC for 6pch and CD for 4pch. The use of Rl-based analysis developed in this 
work on a broader range of patchy systems will be a useful tool to characterize seemingly 
structureless aggregates, and capture early precursors of order. 
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